Load Packages for Mapping

library(sf)
library(dplyr)
library(ggplot2)
library(leaflet)
library(scales)
library(ggmap)

Read in Map Data

ak_regions<-read_sf("shapefiles/ak_regions_simp.shp")
class(ak_regions)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
head(ak_regions)
## Simple feature collection with 6 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.2296 ymin: 51.15702 xmax: 179.8567 ymax: 71.43957
## CRS:            4326
## # A tibble: 6 x 4
##   region_id region      mgmt_area                                       geometry
##       <int> <chr>           <dbl>                             <MULTIPOLYGON [°]>
## 1         1 Aleutian I…         3 (((-171.1345 52.44974, -171.1686 52.41744, -1…
## 2         2 Arctic              4 (((-139.9552 68.70597, -139.9893 68.70516, -1…
## 3         3 Bristol Bay         3 (((-159.8745 58.62778, -159.8654 58.61376, -1…
## 4         4 Chignik             3 (((-155.8282 55.84638, -155.8049 55.86557, -1…
## 5         5 Copper Riv…         2 (((-143.8874 59.93931, -143.9165 59.94034, -1…
## 6         6 Kodiak              3 (((-151.9997 58.83077, -152.0358 58.82714, -1…
st_crs(ak_regions) #Provides projection information
## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["GCS_WGS_1984",
##     DATUM["WGS_1984",
##         SPHEROID["WGS_84",6378137,298.257223563]],
##     PRIMEM["Greenwich",0],
##     UNIT["Degree",0.017453292519943295],
##     AUTHORITY["EPSG","4326"]]

Plotting Data

plot(ak_regions)

#Doesn't look correct, need to transform

Transform Map Projection

ak_regions_3338<- ak_regions %>%
    st_transform(crs = 3338)

plot(ak_regions_3338) #Plot looks much better!

Read in Population Data

pop <- read.csv("shapefiles/alaska_population.csv")

pop_4326<-st_as_sf(pop,
                  coords=c("lng","lat"), 
                  crs= 4326, #set this here according to what the data are 
                  #currently in, not what you want them to be in
                  remove=FALSE)
#head(pop_4326)

#Cannot join pop with ak until they have the same coordinate system

pop_3338<-pop_4326 %>% 
  st_transform(crs=3338)

pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)

plot(pop_joined)

pop_region<-pop_joined %>%
  as.data.frame() %>%
  group_by(region) %>%
  summarise(total_pop=sum(population))

pop_region_3338<-left_join(ak_regions_3338,pop_region,by="region")

plot(pop_region_3338)

pop_mgmt_338 <- pop_region_3338 %>% 
  group_by(mgmt_area) %>% 
  summarize(total_pop = sum(total_pop),do_union=FALSE) #do_union=FALSE maintains the regions delineated
## `summarise()` ungrouping output (override with `.groups` argument)
plot(pop_mgmt_338)

rivers_3338 <- read_sf("shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)
## Coordinate Reference System:
##   No user input
##   wkt:
## PROJCS["Albers",
##     GEOGCS["GCS_GRS 1980(IUGG, 1980)",
##         DATUM["unknown",
##             SPHEROID["GRS80",6378137,298.257222101]],
##         PRIMEM["Greenwich",0],
##         UNIT["Degree",0.017453292519943295]],
##     PROJECTION["Albers_Conic_Equal_Area"],
##     PARAMETER["standard_parallel_1",55],
##     PARAMETER["standard_parallel_2",65],
##     PARAMETER["latitude_of_center",50],
##     PARAMETER["longitude_of_center",-154],
##     PARAMETER["false_easting",0],
##     PARAMETER["false_northing",0],
##     UNIT["Meter",1]]

MAPS!

ggplot()+
  geom_sf(data=pop_region_3338, mapping = aes(fill = total_pop)) +
   geom_sf(data = rivers_3338, aes(size = StrOrder), color = "black") +
  geom_sf(data=pop_3338,mapping=aes(),size=0.5)+
  scale_size(range=c(0.01,0.2),guide=FALSE)+
  theme_bw() +
  labs(fill = "Total Population") +
  scale_fill_continuous(low = "gray90", high =  "purple", labels = comma)

Basemaps

pop_3857 <- pop_3338 %>%
  st_transform(crs = 3857)

# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))
  
  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
  
  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}
bbox <- c(-170, 52, -130, 64)   # This is roughly southern Alaska
ak_map <- get_stamenmap(bbox, zoom = 4)
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)
ggmap(ak_map_3857) + 
  geom_sf(data = pop_3857, aes(color = population), inherit.aes = F) +
  scale_color_continuous(low = "gray90", high =  "purple", labels = comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Leaflet Practice!

epsg3338 <- leaflet::leafletCRS(
  crsClass = "L.Proj.CRS",
  code = "EPSG:3338",
  proj4def =  "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
  resolutions = 2^(16:7))
pop_region_4326 <- pop_region_3338 %>% st_transform(crs = 4326)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = "gray",
                    weight = 1)

pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)

m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
        addPolygons(data = pop_region_4326, 
                    fillColor = ~pal(total_pop),
                    weight = 1,
                    color = "black",
                    fillOpacity = 1,
                    label = ~region) %>% 
        addLegend(position = "bottomleft",
                  pal = pal,
                  values = range(pop_region_4326$total_pop),
                  title = "Total Population")
m